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A  new  approach  to  a  posteriori  error  estimation  is  outlined  which  is  applicable  to  general  h-p  finite 
element  approximations  of  general  classes  of  boundary  value  problems,  I'he  approach  makes  use  ol 
duality  arguments  and  is  based  on  the  element  residual  meihorl  (t:RM)  Important  aspects  of  the 
method  are  that  it  provides  a  systematic  approach  toward  deriving  element  boundary  conditions  for  the 
HRM;  it  leads  to  an  tipper  bound  for  the  global  error  :in  approtsrCo.e  eneigx  iiomi.  and  it  is  valid  for 
non  uniform  and  irreguiar  h-p  meshes.  In  the  present  exposition.  ;i  brief  outline  of  the  theoretical 
foundations  of  the  method  is  given  tvrgcther  with  the  results  of  its  application  tt>  several  represemative 
problems.  These  results  show  that  the  appntach  is  applicable  to  general  linearlv  elliptic  systems, 
including  unsymmetrieal  operators,  and  that  the  method  is  valid  for  broad  chisses  of  linear  and 
non-linear  problems. 


1.  Introduction 

In  a  recent  paper  [1],  we  developed  a  ^icncral  theory  tor  a  posteriori  error  estimation  which 

has  the  following  attributes: 

(1)  it  employs  a  special  variant  of  the  clement  residual  method  [1-.'^]; 

(2)  under  mild  assumptions,  it  produces  estimates  in  convenient  energy  type  norms  which 
may  not  be  directly  associated  with  the  actual  bilinear  form  of  the  problem  under 
consideration; 

(3)  it  employs  a  local  duality  argument  that  leads  to  a  guaranteed  global  upper  bound  to  the 
error  and  which  generalizes  the  duality  method  of  Kelly  [4); 

(4)  it  is  valid  for  symmetric  and  unsymmctric  operators; 

(5)  under  additional  assumptions,  the  approach  can  lead  to  asymploiicafly  exact  error 
estimators; 

(6)  it  is  well  suited  to  irregular  meshes  with  non-uniform  h-p  finite  clement  approximations 
and  functions  independently  of  the  order  p  of  the  local  element  shape  functions: 

(7)  it  employs  a  systematic  scheme  for  flux  balancing  on  clement  boundaries  that  substantially 
increases  the  quality  of  the  local  and  global  cffcctivity  indices. 


(  orrespondeiice  to:  Mark  Ainsworth.  Texas  Institute  for  Computational  Mechanics,  WRW  ,3(b.  The  Universitv 
of  Texas.  Austin.  TX  78712,  U.SA. 

'  r)n  leave  from  Mathematics  Department.  Leicester  University,  t.eiccstcr,  UK. 
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\l  .  1.  I  (  lilcil .  .  \  fUDi  Iilllh'  /(!/  il  jUiMrnntl  triiii  f'.UnttlU'  n! 


The  theory  generali/es  picvitnis  work  on  the  I- RM.  In  [Ktrtieuhtr.  tor  the  special  case  ot  linear 
triangles  { />  =  1).  the  eonjeetiire  inatle  hy  liank  and  Weiscr  (2|  is  conlirmed:  that  a  certain 
variant  i>t  the  I:RM  always  provirles  an  upper  hound  on  the  error. 

Our  purpr)se  in  the  present  paper  is  to  briefly  outline  the  principal  ieatures  oi  the  theorc  in 
connection  with  a  simple  model  elliptic  boundary  value  problem  and  to  locus  on  issues  ol 
implementation.  I'he  robustness  and  generality  ol  tiie  metlu'd  is  ilemonstrated  In  some 
applications  including  elliptic  systetns  and  unsvmmetrica!  problems.  I  he  applicabihtc  ot  the 
method  to  arbitrary  //-/;  meshes  is  a'so  illustrated.  In  partictilar.  it  is  demonstrated  that  the 
method  yields  ver\  good  local  estimtites  both  lor  meshes  with  odd  and  with  even  order  shape 
[unctions  on  neighbrrring  elements,  in  contrast  to  other  technic|ues  proposed  in  recent 
literature. 


2.  Theoretical  foundations 

2. 1 .  Model  problem 

For  clarity,  we  begin  by  considering  a  simple  model  elliptic  boundarv  value  problem  in  two 
dimensions:  Find  u  ~  u[.\.  y)  such  thtit 

+  h  Vii  +  eii  =  f  in  (i  , 

I '  > 

(I  “  =  g  on  1 .  ,  If  =  d  on  / ,,  , 
a/f  ^  " 

w'hcre  12  is  a  connected  Lipschitzian  domain  in  R'  with  boundtiry  ell  -  In  (  1 ),  the 

coctticients  a,  h,  c  and  data  /,  g  arc  assumed  to  be  such  that  n  c.xists,  is  unitjue.  is  cr>ntinuous 
i>n  the  interior  ol  and  depends  coniinuously  on  the  data  in  appropriate  norms.  I  he  weak 
form  of  ( 1 )  is  as  follows:  Find  i<  G  /  such  that 


where 


Bin.  v)  -  L(v)  Vv  E  I  . 
t  =  ( V  6  //  '(fi  ):  yn  =  d  on  I],  j 


and  B  :  f  x  /  — >R.  /_  ;  /  are  the  forms 


Bin.  v)  -■=  I  (r/V/f  ■  Vr  f  vh  Vn  +  env)  dx 


I-iv)  -  fv  dx  +  ge  d.v  . 


2.2.  Parlilioiiinf^ 

We  next  introduce  a  partitioning  of  (2  into  N  ^  .V{.^)  subdomains  (finite  elements)  22^ 
where  fi  =  U  and  construct  on  a  space  1  E  1  of  piecewise  polynomial  functions.  The 
spttce  /  could  include  arbitrary  h-p  finite  element  approximations  of  the  type  discussed  m  |.^1, 
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Following  Standard  finite  element  proeedurcs.  we  supp(.)se  that  a  liinetion  //  G  /  is  obtained 
whieh  represents,  in  some  sense,  an  approximation  of  the  st)iution  u  of  ( 1 ).  Our  goal  is  to  use 
the  available  information  to  ealeulate  an  estimate  of  the  error 

e  ==  II  II 


in  an  appropriate  norm 

With  regard  to  the  partition  we  introduce  the  following  notation;  are  loeali/a- 

tions  t)f  the  forms  in  (4)  and  {.^). 


/^,  (h.  v)  =  (iiVii  Vp  +  vb  Vu  +  CKL’)  djr . 


Bill.  V)  =  L  Bf^lii.  V)  . 
K  I 


/u  djf  +  gu  d,v  . 


\ 


/V 


L{v)  =  S  ■ 


(7 


(H) 


(9) 


K  I 


for  u.  vE  (.  Further  ^  is  a  local  subspacc  of  f.  with 


/  =  0  . 
K  I  ^ 


iO) 


There  now  arises  the  issue  of  the  norm  in  which  we  shall  measure  the  error  c.  For  this 
purpose,  we  introduce  a  symmetric,  positive  definite  bilinear  form  a  \  ‘t  k  'f  —*R. 


(N) 


{i{u.  e)  =  (rtVff  •  Vu  +  cuv)  dx  . 


where  a  and  c  are  constants  which  arc  arbitrary  except  for  the  retniirement  that  the  original 
bilinear  form  B(  - .  ■ )  of  (4)  is  coercive  with  respect  to  the  norm  induced  by  u{  - .  •  ).  1  hat  is. 


where  (3  >0  is  a  constant,  and 

IkL  =  V(tJv.  u) . 


(12) 


(13) 


In  wc  write 

v 

a(u.  u)  X  d.  (i{.  u)  . 


K  1 


(14) 
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a  Ml.  v)  =  (a\H  Vc  +  cuv)  djr .  ( 1?  ) 

Ik’ll'L  =  1’) . 

llt’lL  =  X 

K  I 

We  also  introduce  the  averaged  local  flux 

(«^  •  uW/)(,v)  •  l(  1  -  +  tt^,(.v)uVuj J  .  (18) 

when'  v  is  a  point  on  the  interelement  boundary  ,  n  shared  by  neighboring 

elements,  is  the  unit  vector  exterior  and  normal  to  and  is  a  linear  (unction 

associated  with  edge  1'^,. 

With  this  notation  now  established,  we  consider  the  following  local  problem:  Find 
such  that 

a^{<b^.  w)  =  w)  +  tp  ^  ■  «Vh){.v)  d.v  ,  h’  ^  ^ 

Equation  (19)  characterizes  the  local  problem  providing  the  basis  for  the  error  residual 
method  corresponding  to  the  norm  H-H,.  The  significance  of  (19)  is  embodied  in  the 
following  theorem, 

THEOREM  I.  Let  (fr^,  he  the  solution  of  the  local  problem  (19)  corresponding  to  the  element 
.  Then  the  functions  of  (IS)  can  he  chosen  so  that 


S  \Mf,. 

P  K  1 

where  fi  is  the  constant  appearing  in  (12). 


PROOF.  Sec  (1.6). 

REMARK  1.  In  the  case  in  which  the  bilinear  from  /?(■,-)  is  symmetric  and  positive 
definite,  we  can  take  a{- .  ■  )  =  6(  • ,  •  )■  Then  the  constant  p  =  1  and  the  norm  ||  •  || ,  reduces 
to  the  standard  energy  norm. 

REMARK  2.  The  introduction  of  the  symmetric  form  «( • ,  •  )  is  equivalent  to  symmetrizing 
the  problem  [7], 

REMARK  .1.  The  conditions  on  the  approximate  solution  u  (that  uG  7  )  can  he  weakened 
considerably.  Let  t/y  be  a  standard  degree  1  basis  function  (a  pyramid  function)  associated  with 
a  regular  node  in  Then  we  need  only  require  that  u  E.  C"(/2)n  //  (.7’)n  f  and 


mu.  til)  =  L(tlr)  . 


(21) 
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That  is.  H  need  only  satisfy  an  orthogonality  eondition  vvith  respect  it)  the  lowest  degree  basis 
function. 


REMARK  4.  Let  =diam(i'i^  )  and  suppose  that 
|/?(n,  t')l  ^  ,  llyjl ,  Vn.t'G  f  . 

As  an  additional  assumption,  suppose  that  the  following  condition  holds: 

^  ^  f  o  I  ^  ^  ^ ii ‘'^1 • 

K  .  i 


(22) 


(23) 


where  p  =  p(h.  p)  is  a  constant  possibly  depending  on  the  mesh  parameters  h  and  p.  I'hen  it 
can  be  shown  [1,2)  that 


X  ||d.J|-,,  ^M{l  +  O(/0+Cp}||t>||;  . 


(24) 


K  1 


where  It  =  max^  /i^  .  Then  we  have 


/3||«ir, -5  2  mjIV.sAyiiai;  ,  (25) 

K  ■-  i 

where  M  depends  on  M,  h  and  p.  This  result  establishes  the  equivalence  of  the  global  a 
posteriori  error  estimate  to  the  true  error.  Moreover,  if  p— >0  as  //—»()  then  we  have  M—f  M. 
This  shows  that  the  constants  appearing  in  (25)  are  asymptotically  optimal.  In  the  case  of 
B{  - .  •  )  symmetric  and  positive  definite,  wc  have  asymptotic  exactness  of  the  error  estimator. 


3.  Implementation 

The  actual  computation  of  the  error  estimator  may  be  thought  of  as  consisting  of  two 
distinct  stages: 

(a)  the  calculation  of  the  linear  splitting  function  a^,  used  in  (18)  to  obtain  the  boundary 
conditions  for  the  local  problem  (19), 

(b)  the  (approximation  of)  the  solution  of  the  local  problem  (19). 

The  fundamental  criterion  determining  the  choice  of  flux  splittings  used  in  the  average  (18),  is 
the  following; 


fi,.(«,  1)-L,(1)  + 


•  «W<)(.v)  d.v  . 


(26) 


A  simple  physical  interpretation  of  (26)  is  seen  in  the  special  case  <;/(x)^  l,  ftf-v) --d  and 
c(jir)^0.  For  this  situation  (26)  becomes 


0  =  j\x)  dx  +  I 

JK  JSKfOS. 


g(.s)  d.v  + 


«W<)(,v)  d  .v 


(27) 
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'I'hc  condition  (27)  proclaims  that  the  data  for  die  local  piobicin  (  Ib)  is  in  et|iiilihtium.  or  that 
the  tluxes  have  been  equilibrated. 

Kelly  |4|  used  this  same  criterion  to  determine  Ihix  splittings  m  the  case  of  picccsvisc 
bilinear  linite  element  approximation  of  Poisson's  equation  in  two  (.limensions.  Our  approtich. 
while  rcialed  to  Kelly's,  differs  significantly  in  several  ways.  The  splittings  used  in  our 
algorithm  are  linear  functions,  as  opposed  to  constants:  our  splittings  can  be  obtained  using 
only  local  computtitions  rather  than  iippiying  a  i^lohal  optimization  procedure;  and,  under  mikl 
assumptions,  achieve  the  equilibration  exactly  (subject  to  rounding  error)  and  are  applicable 
to  general  linear  elliptic  systems  of  second  order  partial  differential  equations.  In  addition,  our 
approach  applies  to  general  h-p  linite  element  iipproximations  on  irregular  mesiics.  v\iih 
non-uniform  j>  and  is  valid  for  triangular  elements,  quadrilateral  elements  or  iiulccd  combina¬ 
tions  of  the  two.  Importantly,  our  approach  applies  equally  well  to  one.  two  or  three 
dimensvons.  ) 

3.1.  The  Jln.x  splitting  algorithm 

The  complete  details  involved  in  deriving  the  algorithm  to  he  presented  can  he  found  in  [8], 
Here,  we  restrict  tmrselves  to  the  bare  essentials  necessary  to  implement  the  algorithm. 

A  key  role  in  the  algorithm  is  played  hy  the  degree  one  basis  functions  (that  is.  the  pyramid 
functions  associated  with  the  regular  nodes  in  the  partition).  Denote  the  regular  nodes  by 
A.  B.  .  .  .  and  let  ih  ^  denote  the  pyramid  function  associated  with  node  A  (settled  so  that 
i//,  =  1  iit  node  A). 

The  computations  are  localized  using  the  patches  of  elements  over  which  the  functitins  li/, 

htivc  non-zero  values.  For  ease  of  nottition.  vve  suppose  elements  fi,.  SL . /i.,  constitute 

the  patch  .V  ,  of  elements  on  which  i// ,  does  not  vanish.  Sr)mc  examples  of  possible  ptitches  are 
shown  in  Figs.  l-.l. 

Associated  with  each  patch  .V  ^  is  ti  matrix  7'  ,.  't  he  matrix  depends  only  on  the  topology  of 
the  patch  S  ,  and  not  on  the  geometry.  For  this  reason  we  refer  to  T  ^  as  the  topology  matrix 
for  the  patch.  ,  is  a  square  matrix  of  size  N  x  /V.  where  iV  is  the  number  of  elements  in  the 
patch.  Therefore.  7',  are  typically  sm;ill  matrices  whose  sizes  do  not  increase  as  the  partition  is 
refined.  The  entries  of  /' ,  are  given  by 


r  2  (I  -I  -It 


t  ig.  2.  topology  matrix  for  bouiiciarv  noitc. 


f  ig.  I.  topology  malrix  lor  interior  node  on  regular 
mesh. 
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[ . TcTri 

!  D,  — ! 

/  o, I  o,  ; 

{ . ‘  H 

\  o,  \  a.  j 

I _ L - ‘ 

f'ig.  A,  Topnlugv  matrix  for  interior  nolle  on  1 -irregular  mesit 


1  -1 
t  1 


(Fig.  2), 


2  -1  0  0  0  -1 
- 1  3-1-1  0  0 

0-1  2-1  0  0 
0  -1  -1  4  -1  -1 

0  0  0  -1  2  -1 

- 1  0  0-1-1  3 


(Fig.  3)  . 


'C,.  if/=/. 

(7’^)  :=  .  -1  ,  if  and  arc  neighbours  in  the  patch  ,  (31) 

[o,  otherwise. 

where  C\  is  the  number  of  elements  in  the  patch  which  share  an  edge  with  element  /2,,  Some 
examples  of  topology  matrices  for  various  types  of  patch  in  two  dimensions  arc  shown  in  Fias. 
1-3. 

The  singular  matrix  7,  is  then  modified  by  adding  1  to  every  entry,  thereby  giving  a  now- 
matrix  T  ^  with  entries 

C, +  1,  if/  =  /, 

(7,),^  =  •  0  .  if  42,  and  42^  arc  neighbours  in  the  patch  .  (32) 

1  .  otherwise  . 

where  C,  is  as  before.  It  may  be  shown  [8]  that  7 ,  is  non-singular.  In  fact  it  is  symmetric, 
positive  definite  and  ainscquently  simple  to  invert  numerically  (using,  for  example,  an  LU 
factorization). 

Define  the  /V- vector  b  with  entries 


where 


(*.i)k  =  Sx(“’  */'»)  +  f  .  „  •/4i('>')<«x  •  I  ’  d.s  . 

JrtK 

{ ■  aWt ) , ;  =  .  .v  e  . 


Furthermore,  for  every  interelement  edge  l\j  within  the  patch  define 


where 


■  «Vw[]d.v . 


\n  ■  flT?/]]  ■=  n^.  ■  +  n,  ■  a^u\,  .  s  E  . 


A/.  A'm.wiorih.  J  1 .  Oden.  .4  procednn'  tor  a  ptisiennii  error  estinumon 

Having  calculated  these  variinis  quantities,  a  set  of  consuinis  is  computed  using  the 
procedure  outlined  in  Fig.  4.  The  case  when  .  vanishes  requires  a  little  extra  care,  details 
of  which  can  be  found  in  [8|. 

The  procedure  in  Fig.  4  is  applied  to  every  regular  node  in  the  mesh.  The  actual  tlux 
splitting  used  in  (18)  is  then  given  by 

«a7.('')  =  S  ■  v  e  I'k,  ■  (37) 

T 

Most  of  the  terms  in  this  summation  vanish  due  to  lA,  having  non-zero  values  on  a  small 
number  of  edges.  For  example,  in  the  case  of  regular  meshes,  only  two  terms  in  the 
summation  are  non-zero,  namely  those  corresponding  to  the  two  nodes  forming  the  endpoints 
of  the  edge  l\,  .  i.e.  is  then  the  linear  function  which  interpolates  to  ,  and  „  at  the 
endpoints  of  the  edge.  In  the  case  of  irregular  meshes  the  situation  is  more  complicated  with 
at  most  three  non-zero  terms  appearing  in  the  sum. 

It  can  be  shown  {8j  that,  with  this  choice  of  splitting,  the  condition  (2b)  is  satisfied.  The 
process  described  may  appear  elaborate.  However,  the  computational  work  entailed  is  rather 
small  by  comparison  with  the  cost  of  performing  other  standard  tasks  in  the  finite  element 
method.  In  [8|,  an  operation  count  shows  that  the  process  is  of  optimal  order,  increasing  only 
linearly  with  the  number  of  elements  in  the  partition. 

3.2.  Appro.ximiuion  of  local  problems 

The  approximation  of  the  local  problem  ( 19)  is  performed  by  means  of  a  Galerkin  method 
using  a  particular  set  of  triiil  functions.  Here,  we  shall  describe  the  procedure  we  use  for 
quadrilateral  elements. 

Let  {P„{x)}  denote  the  usual  Legendre  polynomials  on  [-1.  Ij.  It  will  be  necessary  to  be 
able  to  compute  the  values  of  the  polynomials  themselves  and  their  derivatives  efficiently. 
Unfortunately,  in  many  textbooks  it  is  suggested  that  these  quantities  be  calculated  directly 
from  the  expansion  in  terms  of  powers  of  .v.  This  approach  is  not  only  unnecessarily  expensive 


for  each  regular  node  A  do 
begin 

calculate  ; 

calculate  an  LU  factorization 
for  every  element  K  E  do 
begin 
calculate 
calculate  Pk/'.a 
end; 
solve 

for  every  interelement  edge  I].,  in  the  patch  do 
begin  . 

i  ,A~  ^L,A  ^  /  PkI,,A 

end 
end ; 


Fie.  4.  PscikIo-coiIc  i)f  llu.x  splitting  algorithm. 
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but  leads  ti)  catastrophic  build  up  ol  rouful-tdl  errors.  'HK-rc!orc,  we  use  the  three  term 
recurrence  relations,  f'hus  ,  lor  st)me  lixeil  value  ot  .v.  vve  calculate  /’„(  x  )  aiui  O  as  tollows: 


/>,(.v)-,v.  t.'S) 

.  i(-v)  -  l\  ,(  0  .  •  I  . 

and 

SyT  4-  1  A  i-  1 

.  ,(.v)  =  .v/’;(,v)  -  -  Y  r:  ,(.v) .  a  •  i . 

To  calculate  the  values  of  /'.,(.v)  and  P'„(x)  in  this  way  requires  oitly  order  r,  operations. 
Mrmeover,  as  a  byproduct,  one  ;ilso  obtains  the  values  ot  all  the  lower  order  polynomials  and 
their  derivatives  at  no  extra  expense. 

The  values  of  l\^{x) . f\,ix)  and  P^ix) . /V,{ ' )  are  then  used  to  compute  the 

functions  ,v„(v) . A',,!-')  ‘‘'"'‘J  A'h(') . v',(')  given  by 


A-hI-v)  =  \  \  /*,{-v)  .  A,(-t)  ■"  \  i  .  yd.v)  =  ^  \  --  /’.(.V)  . 

_ _  '  ■  “  (40) 

A'xl  O  -  Y{k^)  ^  ^  . "  ■ 

and 

V =  <*  •  4  \  i  /'„(•'■)  .  A':(-V)  =  \  PAx)  . 

"  ■  <-n) 

a;(.v)  -  Y  -L--  /Y  ,(.v)  .  A  -  3 . n  . 

Let  fi  =  [  -  1 .  1 1  X  I  ~  1 .  1 1  denote  the  usual  reference  element.  The  trial  functions  will  be 
detined  on  fi.  For  simplicity,  assume  that  the  linite  element  approximation  u  is  a  complete 
VX'lynomial  of  degree  p  on  the  element  under  consideration.  Let  q  -O  be  an  integer;  then  we 
delinc  the  space  L '  ''(/))  by 


T'’ ''(/i)  ^  {  A,(^)Ta(^)  I  of /.  A  .>/;}  .  (42) 

It  is  seen  that  dirn(T'' '7  =  (/->  + r/ +  1  )"-(/■'+  1  )“■  The  index  q  controls  the  number  of 
increments  in  the  polynomial  degree  of  the  space  and  may  be  used  to  increase  the  dimension 
of  the  space. 

I'he  local  apprr)ximation  space  is  then  taken  to  he  ^  D  V'”'' ).  The  discretized 

version  of  the  local  problem  (19)  is:  Find  G  ‘f  ^  such  that 


(Y  .  ^x)  =  T^(n')  -  +  9  'x{n^  ■  aVu){s)  d.v  Viv  G  / 


■>I.K 


(4.7 


Owing  to  the  above  assumptions  and  the  construction  of  the  local  space,  this  purblem  always 
has  a  unique  solution. 
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4.  Numerical  examples 

In  this  section  wc  present  three  examples  to  illustrate  the  perlormance  oi  the  algorithm. 
In  order  to  compare  the  estimated  error  with  tire  true  error  it  is  necessary  to  accurateis 
calculate  the  true  error  over  each  element.  In  all  our  examples,  the  true  error  is  computed 
using  an  algorithm  which  approximates  the  integral  using  lirsl  a  single  subdomain  and  then 
subdivides  the  region  rd  integration  into  tour  subdomains  and  approximates  the  integral  o\er 
the  lour  subregions.  IWo  estimates  are  thus  obtained  tor  the  true  value  of  the  integral.  If  the 
difference  in  these  two  estimates  is  less  than  \'  (  relative  error,  then  the  approximation  is 
accepted.  Otherwise,  the  element  is  further  subdivided  into  lb  regions  and  so  on  untii 
agreement  between  consecutive  approximations  is  obtained  to  less  than  lO  relatixe  errem. 

4.1.  .Svnunerric  eUijUi^.  i^pcrator  u  ith  .sniooih  .solution 
rhe  problem  we  consider  is;  F  ind  ti  such  that 

Sit  u  ~  i)  in  <}  .  (44 ) 

subject  to  tlie  biHindary  exmditions 

n(().  y )  “■  11  sin  iry  .  (I  v  1  /2  . 

(45) 

n(.v.d)  -  n(.v.  1 /2)-(l  .  ()<.v<l/2.  ^--0.  .v  -  I /2.  0  •  v  ■  I /2  . 

<tn 

where  <r  ■■■■■  2Tr.  The  geometry  of  the  domain  and  the  boundary  couvlilions  applietl  are  shown 
in  I  ig,  5. 

I'he  true  solution  is  given  by 

ii{.\\  V)  --  (cxp|(.v  -  1  )V  1  -f-  ir'l  exp(  \\  1  i-  rr'  ))  sin  trv  .  (46) 

In  this  case,  wc  have 


HUi.v)  --  j  (Vn-Vr  t  nrd  dx  (47) 

and  wc  choose  a(n,v)-  B(ii.v).  rheorem  1  predicts  that  we  will  obtain  guaranteed  upper 
bounds  on  the  true  error  measured  in  the  energy  norm  defined  by  ■  .  •  )•  provided  that 
equilibration  of  the  tiuxes  is  achieved. 

Ihe  problem  is  solved  using  uniform  meshes  of  quadrilateral  elements  with  uniform 
polynomial  degree,  The  local  problems  are  a|iproximaled  using  an  increment  r/  -  2  in  the  local 
space.  That  is.  a  degree  p  finite  element  apprrrximation  is  ;maly/ed  using  the  space  V  '■  .  In 
each  case,  the  equilibration  procedure  is  able  to  reduce  the  lack  of  equilibration  to  the  level  of 
round  off  error  on  virtiuilly  every  element  in  the  partition. 

I'able  1  contains  the  results  obtained  for  finite  element  approximations  of  degree  1-4  on 
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Number  oi  elements 

lb 

!  2S 

1 

t  .01  :‘)IW 

1  .Utlt.s.sS 

i  nniipir 

■> 

1  .Olio’s?: 

1  .iiii:iti'i 

1  lllllltbll 
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I .  bill  II  IPS 

i  .UIIOIll'i 

4 

1  III !( 0.47 

i  .iiiiiii:i) 

1  III  llill.'l) 

meshes  etmtaining  16.  64  and  12<S  elemenls.  The  quantity  shown  in  the  table  is  the  cl'firiivity 
index  (the  ratio  of  the  estimated  error  to  the  true  error).  Theorem  I  predicts  that  the 
efieclivitv  index  be  greater  than  unitv.  This  prediction  is  borne  out  bv  the  results  shown  in 
labie  1.  " 

4.J.  (  racked  panel  prohle/n 

Consider  the  problem;  Find  a  such  th;tt 

hi  0  in  , 

subject  to  the  boundary  conditions 

ii(r.  rr)  =  0  ,  0  <  r  <  1  .  an! i)n  -  d  ,  0  <;  r  <  1  .  =  0  .  (49) 

with  u(r.H)-r'  cos  'ff  on  the  remaining  portion  of  the  boundary.  The  geometry  of  the 
domain  <l  is  shown  in  Fig.  6. 

The  true  solution  is  given  by 


I 

{4S) 


u(r.  ff )  =  r'  '  cos  [O  .  ( 5(1) 

The  problem  is  the  analogue  of  a  cracked  panel  problem  in  linear  elasticity,  vv'th  a  singularity 


u  =  t’'2cos4!? 


^ - u  =  0 


.  . V 

rin 


f  ie  (v.  (ieoincirv  aiut  bouiKlarv  coiKiiUons  lor  crae- 
ked  panel  problem. 
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at  the  orimii.  In  this  casi',  we  have 

/^(  (/,  r  )  \u  •  Xv  lij  (51) 

Ji! 

and  we  take  ii(u.  t  )  r).  I  heorem  I  av;ain  prediets  that  we  will  tihtam  upper  hounds  on 

(he  true  error.  proMdetl  that  we  eiimlihrale  the  lUixes. 

In  analy/inti  this  prohlein,  some  eare  must  he  taken  with  the  approximation  ot  the  loeal 
prohiems,  1  heoretn  1  assumes  that  the  toeal  problems  are  sotxed  exaetlv  ,  which  is  not  tile  case 
in  practice.  Ihererore,  lor  the  purposes  td  illustrating  the  theorx.  a  sequence  oi  approxi¬ 
mations  to  the  true  solution  of  the  ioc.il  problem  is  ohtainei)  h\  incrementing  </.  That  is  to  sa\ . 
we  compute  a  sequence  of  ac»proximations  using  the  spaces 

r-"  '  c:  v"  '  c  r"  ■■  c--- 

until  the  slifference  m  the  norm  of  the  .ipproximation  is  suflicientK  small. 

One  other  feature  of  this  particular  example  is  the  diflicultv  in  estimating  the  error  m 
elements  adiaeent  to  (he  singular  point.  Our  theory  makes  no  promises  concerning  the 
effecti\it\  of  the  estimator  in  a  single  element.  However,  we  present  lesults  slniwing  the 
estimated  and  the  true  error  m  the  elements  adjacent  to  the  singularity  fi,  ami  (see  l  ig.  7). 

lahles  2-4  contain  results  of  estimating  the  ern>r  in  the  approximation  ohtameil  using  a 
unilorm  mesh  of  .'2  quadrilateral  elements  with  uniform  polynomial  degree  1.  l  ive  increments 
m  the  local  approximation  space  are  neeiied  before  agreement  is  obtained,  l-or  purposes  of 
comparison,  we  also  give  the  results  obtained  when  no  equilibration  or  balancing  of  duxes  is 
pertormei.1  (instead,  a  simple  averaging  is  applied). 

Tables  2  and  .5  contain  the  error  estimates  in  the  elements  <idjacent  to  the  singularity  li,  and 
Slf..  The  estimate  ot  the  error  obtained  when  equilibration  of  duxes  is  performeit  is  superior  to 
the  estimate  oldained  using  a  simple  averaging.  Table  4  contains  the  global  estimates  of  the 
error.  It  is  seen  that  it  is  necessary  to  approxin..  te  the  local  problems  accurately  if  one  is  to 
obtain  the  iqiper  bound  proclaimetl  by  Iheorem  I. 
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Table  4 

Lffeet  of  solving  local  problem  with 
cracked  panel  ( /i  ■  I .  .'2  elemenis) 

iiKTcaMni:  «acui.»c\ 

on  tire  estim.iti. 

s  of  global  error  lor 

I'.stimated  global  error 

(ilobal  e 

flectivitv  indev 
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balancing 

balancing 
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Table  .s 

fdlect  of  solving  local  problem  with 
element  il,  for  crackeii  panel  ( /»  - 

iiKToasing  accuracy 

1 .  12N  elements) 

on  the  estimates  ol  local  error  m 

listimaled  local  error 

I  oc:il  effeetiv  ity 

iiulev 

Number 
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With 
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increments 

balancing 

balancing 

balancing 

balancing 
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Tables  5-7  show  the  ccirresponding  results  obtained  when  the  mesh  is  refined  uniformly  to 
!2H  elements  of  degree  one.  The  results  obtained  are  similar  to  the  ease  of  .^2  elements.  Tables 
X-IO  contain  the  results  obtained  when  the  degree  of  the  elements  is  increased  uniformly  t<i 
degree  2  on  M  elements.  Once  again,  the  results  show  the  superiority  of  the  estimate  obtained 
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0,8027b4(  1  ) 

0.8(I2764(  1 ) 
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using  equilibrated  fluxes.  In  eaeh  case,  the  result  of  Theorem  I  is  verified,  although  it  is 
necessary  to  solve  the  local  prt)blems  very  accurately. 

4.3.  Umymmetric  elliptic  sy.stem 

As  a  final  example,  we  consider  the  unsymmetric  elliptic  system  with  non-constant 
convection  given  by;  Find  sueh  that 
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I'ahle 

mtcct  of  solving  local  problem  with  increasing  accuraev  on  the  estimates  ol  local  eiror  in 
clement  for  cracked  panel  (  p  2.  32  elements) 


Number 

of 

increments 

I-stmiated  local  error 

Local  ellectivitv 

lUilcx 

With 

balancing 

Without 

balancing 

With 

balancing 

Without 

bal.incmg 

1 

()..S.v2473( 

i) 

(».b(i.S535t 

1) 

(l.bSbilOd 

(I..S'7427 

T 

O.hl.^bSUf 

1) 

(l.7|705(v( 

1  ) 

(1  7U(Krv3 

().U23,S(I3 

3 

il,fi2S18S)( 

1) 

(I,84,S‘I‘X»( 

1) 

ll.S()'t3i3 

1,(W3777 

4 

().(v44()7y( 

1) 

l),,S7b.Sb3( 

1) 

l>.82'»7S.s 

1. 12b33u 

5 

().64.S2‘W( 

1) 

(),Ul2l.4n( 

1) 

tl  S3.S222 

1  173122 

I'rue  value 

().77fi2()0( 

1) 

(),77b2(M)( 

1) 

i  .OIKKHHI 

I.IMHKIIM) 

Table  lit 

Effect  of  solving  local  problem  with  increasing  accuracy  on  the  estimates  ol  global  error  lor  crackcvl  panel 
(  p  -  2.  32  elements) 


Number 

of 

increments 

tistiniiited  global  error 

(jlobal  eOectiviiy 

iiulcx 

With 

btilaneing 

Without 

balancing 

With 

bahmemg 

Without 

balancing 

1 

(l.94.S()T3(  •  1 ) 

0,79(H)96(  1 ) 

O.S3S555 

0.6v).s,si0 

T 

0. 1(17495(0) 

0.9U473b(  - 1 ) 

0.950753 

0. SI  >0207 

3 

0.11 15.S0(U) 

0,9,S()423(  1 ) 

0.9SWil.S 

O..S6714S 

4 

0.113117(0) 

0.I(H)715(0) 

I.0(KI47,S 

(I..S9(I7S7 

5 

0.114017(0) 

0.104075(0) 

1 .00,S43S 

0.920,305 

True  value 

0.11.^063(0) 

0.11306,3(0) 

I.OOO(KH) 

L0<MMKNI 

-  f  Att, 

flu, 
+  A-  - ^ 

flu, 

—  V  - 

+  .vya. 

-  a,  =0 

I 

(IX 

fly 

-  F  Aa, 

<1  u , 

T  A  - 

flu. 

~  3  ’  Ty 

-  Ava, 

4-  a,  =  0 

rTV 

'  rTv 

where  f  =  1  /lOO,  subject  to  the  boundary  conditions 

=  cxp((.v'  -  >’■  --  1 )  /f't  ■  "  -W  cxp[(-v'  -  y"  -  1 )  /i-  j  on  / 

and 

dU.  ,  , 

F  - —  =  fn  •  V  exp  (t'  -  y-  -  1 )  /f  , 

()n 

f  =  t  n  Vat  exp((.v‘  -  y'  -  1 )  /f  j  on  l\  , 

where  /i.  1],  and  arc  shown  in  Fig.  8. 

The  true  solution  to  this  problem  is  given  by 

a  I  =  exp((A  ’  -  y’  1 )  /f'|  .  a.  -  -y  cxp({A  "  --  y 


(52) 


(53) 

(54) 


1)/H- 


(55) 
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ilb  (1.5) 


<i.o)  f'N  <’-o> 

F'ig.  S.  (iconictry  ;i!u)  hoiituiary  coiuliiion'.  tor  iit)\yniinciii(.'  cllipiu  s\sictit 

The  main  feature  of  the  solution  is  the  presence  of  a  strong  boundary  layer  effect  along  the 
right-hand  wall  of  the  dt^main  caused  by  the  non-constant  coincction  dominated  (lov\  field 
b  =  (.V,  -  y). 

In  this  case,  the  bilinear  form  B(  ■  ,  •  )  is  unsymmetric.  We  chortse  the  bilinear  ft>rm  </(  •  .  •  ) 
to  be 

a{it.  v)  =  j  ^  ( J-'V/Zi  -  Vc,  +  cVn .  ■  Vc,  -l-  +  n,i’, )  dar  .  ( 5b) 

The  theory  presented  in  |6]  shows  that  the  error  estimator  bounds  the  true  error  measured  in 
the  symmetrized  norm.  For  the  purposes  of  illustration,  in  this  e.xampie  we  compute  the  true 
error  in  the  symmetrized  norm  explicitly.  It  is  this  quantity  which  is  Uihclied  as  the  true  error 
in  Table  11. 

The  presence  of  the  boundary  layer  indicates  that  an  adaptive  finite  clement  analysis  based 
on  refining  the  mesh  and  enriching  the  degree  of  the  approximation  is  suitable.  Fhe  sequence 
of  meshes  generated  during  the  analysis  is  shown  in  Figs.  9.  12.  15.  1<S  and  20.  Fhe  meshes  are 
not  only  irregular  but  contain  elements  of  differing  polynomial  degree.  The  final  mesh 
contains  elements  of  degree  six  near  the  boundary  layer.  Nevertheless,  the  behaviour  t>f  the 
error  estimator  remains  highly  satisfactory  as  shown  by  the  results  in  Fable  11. 

Onc  source  of  concern  when  estimating  errors  for  this  type  of  problem  is  that  the 
distribution  of  the  estimated  error  will  not  agree  with  the  distribution  of  the  true  error  owing 
to  the  convective  effect.  Fherefore,  in  Figs.  9-23,  we  present  plots  showing  the  distribution  of 
the  true  and  estimated  errors.  It  is  observed  that  the  distribution  of  the  estimated  error  closely 
reflects  the  actual  error  distribution. 
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Behaviour  of  error  estimaU'rs  lor  utwymmelrie  elliptie  system 


Mesh 

number 
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of 

freedom 
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error 
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global  error 

(ilobal  cffectivUs  index 
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balancing 

W  ithoiii 
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1 

2.S 
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T 

51 
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A 

111 

0.56.3W2(  1 ) 
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4 

IfO 
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j|.ri49044(  -  2) 

I.IHII2.35 

1  003400 

5 

340 

0. 26.3:3451  2) 
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1.0407,30 
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MESH  1 


D.O.F=  25 


f-'ig.  Adaptive  analysis  of  unsymmctric  elliptic  system.  Mesh  1. 


MESH  1 


MIN=  0.570E-0A 
MAX=0.300.A809 
ERROr’=0.40fll28'^ 
D.O.F=  25 


ig.  it).  Adaptive  analysis  nf  unsymmelrie  elliptic  system.  Distrilnitiim  of  true  error  on  Mesh  1. 


<■«» 


M.  Ainwsinili ,  J  1  Oden.  ,1  prtKcdim'  tin  u  poMcrum  cunt  e^ltiiuiiuin 


MESH  1 


(1 


0.25  0..?:5 


M1N'=  0,235E-03 
MAX=0.3045585 
ERROR=0,420103 
D.O.F=  25 


Fie,  II,  Ai.iaptivc  analysis  dl  unsymmctric  elliptic  system.  Distrihntion  of  estimated  error  on  Mesh  I. 
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D.0.F=51 


Fig.  12.  Adaptive  analysis  of  unsymmctric  elliptic  system.  Mesh  2 
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0  0,03  0.06  0.09  0,12 


MIN=  0,209E-03 
MAX=0, 112570! 
ERROR=0, 145063 1 
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F-ig,  1.3.  Adaptive  analysis  of  unsynimelric  cliiplic  system.  Distribution  of  true  error  on  Mesh  2. 
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M1N=  0.108E-03 
MAX=0. 1 129691 
ERROR=0.1447‘ 
D.O.F=51 


F  ig.  14.  Adaptive  analysis  of  unsymmetrie  elliptic  system  Distribution  of  estimated  error  on  Mesh  2. 
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ERROR=0.056399 
D.0.F=  1 1 1 


lA.  Adaptive  analysis  of  unsymnietric  elliptic  system.  Distribution  of  true  error  on  Mesh  3. 
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- .  ERROR=0.05640:: 
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Fig.  17.  Adaptive  analysis  of  unsymmetric  elliptic  system.  Distribution  of  estimaieti  error  on  Mesh 
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f'ig.  HJ.  Adaptive  analysis  of  imsymmctric  elliptic  system.  Distribution  of  true  error  on  .Mesh  4, 
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ERROR=0.0064764 
D.O.F=  165 


f  ig.  20.  Adaptive  analysis  of  unsymmetric  elliptic  system.  Distribution  of  estimated  error  on  Mesh  4, 
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MESH  5 


D.O.F=  340 


Fig.  2t.  Adaptive  analysis  of  unsymmctric  elliptic  system.  Mesh  .s. 
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MTN=  0.633E-05 
MAX=0.0018429 
ERROR=0.002633 
D.O,F=  340 


Fig.  22.  Adaptive  analysis  of  unsymmctric  elliptic  system.  Distribution  of  true  error  oti  Mesh 
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